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Abstract 

An evolution of radiant shock wave front is considered in the framework of a recently pre- 
sented method to study self-gravitating relativistic spheres, whose rationale becomes intelligible 
and finds full justification within the context of a suitable definition of the post-quasistatic ap- 
proximation. The spherical matter configuration is divided into two regions by the shock and 
each side of the interface having a different equation of state and anisotropic phase. In order to 
simulate dissipation effects due to the transfer of photons and/or neutrinos within the matter 
configuration, we introduce the flux factor, the variable Eddington factor and a closure relation 
between them. As we expected the strength of the shock increases the speed of the fluid to 
relativistic values and for some critical ones is larger than light speed. In addition, we find that 
energy conditions are very sensible to the anisotropy, specially the strong one. As a special 
feature of the model , we find that the contribution of the matter and radiation to the radial 
pressure are the same order of magnitude as in the mant as in the core, moreover, in the core 
radiation pressure is larger than matter pressure. 



1 Introduction 

During the implosion that forms a relativistic compact object, nearly all of its gravitational binding 
energy ((GM 2 )/ R ~ 5 x 10 53 ergs ~ 0.2Mc 2 ) is stored as internal energy of a proto-neutron star and 
driven by neutrino diffusion. Roughly speaking, there is a consensus that this collapsing scenario 
requires, three main "ingredients" pQ: 
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1. a copious emission of radiation, been a consequence of the microphysics of the system, tends 
to abandon the system, but the absorption and the scattering in the medium hinder it to 
escape freely. 

2. phase transitions that can induce local anisotropic pressures (i.e. P r ^ P±)- An increasing 
amount of theoretical evidence strongly suggests that, for certain density ranges, a variety of 
very interesting physical phenomena may take place giving rise to local anisotropy (see [H [3] 
and [U EJ El QH CCD EE1 [[3] for more recent studies ) . 

3. the formation and propagation of a surface of discontinuity: a shock wave, a detonation or 
combustion front (deflagration or slow combustion) having width very small compared to the 
size of the system. 

Although now exist several independent numerical codes which provide accurate modeling of grav- 
itational collapse in full General Relativity (see [T3] for a good review on this subject and/or visit 
some these links concerning codes for simulations 15 ), none of these codes provide all the above 
mentioned "ingredient" . 

In order to explore the influence that the scheme of radiation and the local anisotropy exert 
on the propagation of a surface of discontinuity, we shall follow a "seminumerical" approach which 
starting from a known interior (analytical) static spherically symmetric ( considered as "seed") 
solution to the Tolman-Oppenheimer-Volkov equation. It can be considered as a continuation of 
previous studies 25, 26, 27J. This scheme transforms the Einstein partial differential equations into 
a system of ordinary differential equations for quantities evaluated at the surfaces (boundary & 
shock) . It is an extension of the so called HJR method [TB] , which has been successfully applied 
to a variety of astrophysical scenarios (see [24] and references therein) and which has been recently 
revisited [LZl CH [19] . 

When considering the interaction between radiation and ultradense matter, we consider two 
quantities: the flux factor, f — T / pR, and the variable Eddington factor, \ = Pr/ Pr, and a closure 
relation between them, i.e., x — x(f) ( see HH1 HD HH H3] ), where T, pn and Pr are the radiation 
flux, and the contribution of the radiation to the energy density and radial pressure. 

Here we shall show that the strength of the shock appears to be a very significant feature 
concerning the evolution of the distribution, because it drastically increases the fluid velocity behind 
the shock. We have also found that the energy conditions (specially the strong one) are very sensitive 
to the variation of the local anisotropy, and hydrodynamic & radiation pressures emerge of the same 
order of magnitude and even, at the core, radiation pressure is larger than that of matter counter 
part. 

This paper is organized as follows: in section [2] the energy- momentum tensor and field equations 
for a non-static spherically symmetric fluid with matter and radiation are established. The Post- 
Quasistatic approximation is considered in section [3] Section [5] is devoted to describe the Taub 
junction conditions across the surface of discontinuity. Finally, we end this work with discussing 
some the results of the numerical simulations. 

2 Matter, radiation and field equations 

As it can be appreciated in Figure [TJ the collapsing configuration is described by three regions: a 
core (J) as the more inner region, the mantle (II) in the middle, and the outer space labeled by 
(III). 
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Figure 1: Regions of the matter distribution. The 
core is labeled by I is the more compact one; the 
mantle (77) is the region in the middle and the 
exterior space-time (III). 



Core and mantle are separated by mean of a shock (discontinuity hypersurface) and a boundary 
surface (hypersurface with velocity equal to the fluid in these point) is define the matter configu- 
ration. Both the shock and the boundary surface are time-like hypersurfaces and there only exists 
unpolarized radiation emerging which is modeled by the radiative Vaidya exterior metric. 

The line elements at each region are 

dsj u = e v dt 2 - e x dr 2 - r 2 (d9 2 + sin 2 (9d</> 2 ), (1) 



and 



ds 2 ln = ( 1 - ) du 2 + 2dudR - R 2 (d6 2 + sin 2 6d< r ) . i 2 ) 



where v and A are functions of t and r, u is the retarded time and R is a null coordinate (gnu = 0). 

For a comoving observer with radial velocity to relative to local Minkowskian coordinates, the 
physical content of the fluid is represented by the radiation flux T in the radial direction; the energy 
density p and P ri and P± as the radial and tangential pressures, respectively. The bar stands for 
the total (hydrodynamic + radiation) contribution. Thus, the energy-momentum tensor for each 
region can be written as, 

(T a p)iji = (p + P±)u a u ~ P±g af 3 + (P - P±)XaXp + 2F( a Uf3) 

= (3) 



where 



p = p + pr, P = P r + P R , P ± =P ± + (P ± ) R , (P ± ) r = PR - Pr , (4) 
e^ 2 S° a - Loe x ' 2 5i ~uje^ 2 S Q a + e x / 2 8j 

u a — r 7t ? Xa — r 7r ; ot ->Xu i (5) 

with the subscript R denoting the contribution of the radiation. 

For the above energy-momentum tensor, Einstein field equations G a p = — 87rT"g become 

P + u 2 P + 2ujT\ _ 1 i _ x ( 1 A' 
l^J 2 J ' 



r 2 \ r 2 



«" A 3-T • (6) 
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8ttPj_ 



2A + A(A-;>) 



2 y" + + - (i/ - A') 



uj{p + P) + (1 + J 1 )? 



(8) 
(9) 



1-CJ 2 

where dots and primes are denoting time and radial derivatives, respectively. 

3 The Post— Quasistatic Approximation (PQA) 

The quasistatic regime is "the next step coming out from hydrostatic equilibrium" . In this regime 
the distribution changes slowly in a typical time scale which is very long compared with the char- 
acteristic scale within the sphere reacts to a perturbation. Mathematically, it can be stated as [18] 
O{uo 2 ) = A = v = \v = A 2 = v 2 = . 



Now, we define the "efective variables" 
, 5 + Plu 2 + 2ojT 



1-lu 2 



and 



P 



P + puj 2 + luiT 
~ ~ 2 ' 



(10) 



which satisfy the same set of equations as their corresponding physical variables in the quasistatic 
case. Thus, Herrera et al. [H] define the post-quasistatic regime as that corresponding to a 
system out of equilibrium (or quasiequilibrium) but whose effective variables share the same radial 
dependence as the physical variables in the state of equilibrium (or quasi-equilibrium) , i.e., is the 
closest possible situation to a quasistatic evolution. 

In this effective variables, the field equations ©-© can be written as 



m 



SttP, = 



4:Trr 2 pdr , 

2(47rr 3 P + m) 
r(r — 2m) 



dr . 



2A + A(A-i>) 



2 y" + [ v > + ±)( v >- V) 



47rr 2 e 2 



[w(p + P) + (1 - lu 2 )T] 



l+UJ 2 

where the mass function is defined by e~ A = 1 — 2f "|*' r ^ . 

4 Collapsing Model with Shock wave in the PQA 
4.1 The regions and the equations of state 



(11) 
(12) 

(13) 
(14) 



Following Herrera and collaborators [251 El] , we consider the core equation of state inspired by the 
anisotropic Schwarzschild-like interior solution, i.e. [2 [29] 



Pi = /(*) 



and 



p _„ f3(l-8/37rr 2 p 7 )^/ 2 fc(t)-l 
1 ~ 91 1 3- (1-8/3^^)^/2 k(t) 



(15) 
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where = 1 — 2/i/ is the anisotropic parameter with hi = corresponding to an isotropic model. 
For the mantle we have an anisotropic Tolman Vl-like solution, i.e. 



Pn 



3g(t) 



and 



Pn = 



Pn 



v / 4=367 



1 - 9D(t)r 

i - Dity^u 



(16) 



again, hn (or £jj = 1 — 2hjj) is the parameter messuring the anisotropy of this region and the 
time-dependent functions f{t),k(t),g(t) and D(t), are obtained from boundary conditions. 

4.2 Junction Conditions 

The matching across the boundary surface to the Vaidya metric is obtained though the Darmois- 
Lichnerowicz junction conditions [3D] , while across the shock front via relativistic Rankine-Hugoniot 
conditions [31] . As it was pointed out by Bonnor and Vickers some years ago |30j . Darmois- 
Lichnerowicz conditions are not equivalent to the O'Brien and Synge ones. In the present case, 
O'Brien and Synge conditions lead to a contact discontinuity (not a shock) forcing the fluid velocity 
to be continuos across the surface r = c(t). Thus denoting the boundary surface by S = r—a(t) = 0, 
the matching conditions lead 



a(u) 



and 



Pa(t) = Fa(t) ■ 



(17) 



Now, if the shock front is described by T = r — c(t) = 0, the continuity of the second fundamen- 
tal form implies the continuity of the mass flux across the shock, and the continuity of the first 
fundamental form and the Rankine-Hugoniot conditions, lead to 



s(t) - Wc(t) - , 



( „_ A)/2 f u(p + P) + {l + ^)T \ _ 



cp 



= 



c(t) 



(18) 



and 



(A_„)/a / ^(P + P) + (1 + ^ 2 )^ 



1 - w 2 



P 



0, 



c(t) 



; (t) 



with [X] 

behind. Now, by using the field equation (fl"4|) in equations (|18p and l|19p we obtain 



(19) 



c(t) + 



X\ r=c u)-i where + stands for ahead from the front shock and — for 



[m + 4irr 2 cp} c(t ) = 



and 



•47rr 2 Pl, 



i(t) 



0. 



(20) 



Expanding the mass function in a Taylor series about the shock front, we can get a relation between 
time and radial derivatives and it implies the jump of the time-derivative of the mass across the 
shock, i.e. 

r )\c(t) = rh(t, c) - cm'(t, r)\ c(t) => [m] c(t) = [-4?rr 2 'cp) c{t) , (21) 
and by using ([2"0|) we obtain 



[P]c(t 



(22) 
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4.3 Metric Functions and Effective Variables 

Let us to introduce the dimensionless variables 



M 



m„ 



A 



F = 1 - 



2M 



n = 



and t 



m„(0)' m a (0)' A m a (0) 

where m a (0) is the total initial mass of the distribution, and the physical variables 

M 

E = Ana 2 ^ , L = —M , and E = -—, 

F 



(23) 



(24) 



with E the luminosity for a (non)comoving observer and L the luminosity at infinity. 

All physical variables can be obtained as functions of the above surface variables A, F, and f2. 
Thus, from the boundary conditions (|17p we get 



Mr 

m n = —: E=(l + n)E, 



Pli 



F. 



1-F 
8nr 2 



P, 



II 



l-F 

247T7- 2 



ip - 9x(r/a)^ 4 ~ 3 «" 



F 



iP-xir/a)^^ 
ip-X 



2/V4-3en ' 



(25) 
(26) 

(27) 



where 



V> = 3(3 + fi)(l- F)6(l + fi).E, x= (l + 3n)(l-F)-6(l + n)J5. 
Now, by introducing the shock "force" parameter iV (P/) c = N(Pji) c , we obtain 



Mc, . ,» . 3M 

mi = — (r/G) ' W = 4^I' 



x 2Mc ( 



3u^/ 2 /c(t) - 1 



Ar 



-(r/cY 



3 - u^/ 2 k(t) 
ffu*(3-fcu«'/ 2 ) 8 /«', 



(28) 
(29) 



where 



k = 



u = 1 



24ttc 2 / 5 / /3 + 3JVa(l - F) 
F^/ 2 [12ttc 2 pi(3 + Na(l - F~j] 

1 3(1 — F) 



# = 



/3 

S(F-l) 



iirr 2 pi 



(3-fcF«'/2)8/«, 
a = i/, - 9x(c/A) V4 - 3 «" , /? = V - x(cM) V4_3? " 



3 ' 2 167rc 2 /5/ 

Again, all effective variables and metric functions depend on t through the surface variables 
A, F, n and L. 



G 



4.4 The Surface Equations 



In order to find the evolution of the surface variables, we have to integrate a system of ordinary 
differential equations on A, F, Q and L, i.e. the surface equations. The first surface equation can 
be obtained from the relation between the coordinate velocity and the velocity of the comoving 
observer r = e~z~u , evaluated on £. The second equation emerges from time derivative of F and 
by using the definition of the luminosity L 



A = FCl. 



F 



(1 - F)Fn + 2L 
A 



(30) 



where we have used the junction condition on A and v (fT7| and the dimensionless variables ((23)) . 
The differential equation for f2 is obtained evaluating the conservation law T" Q = 0, on S, i.c 



P 



(p + P)(47rr 3 P + m) _ 2(P± - P) 
r(r — 2m) r 



Anr{r — 2m) 



3m 2 



2m 



mv . 



but we do not write it for its cumbersome form. 

Because this system is overdeterminate, we should introduce one of the "surface" functions. 
Since the only observable quantity entering a "real" gravitational collapse is the luminosity, it 
seems reasonable to provide such a profile as an input for our modeling. Therefore, we select the 
luminosity profile to be a Gaussian pulse centered at t = to 



■ AM rad 
M = L = ;=- exp 



/2tt 



1 ft — to 
2 



(32) 



where s is the width of the pulse and AM ra d is the total mass lost in the process. 



4.5 The radiation hydrodynamic environment 

Conscious of the difficulties to cope with dissipation due to the emission of photons and/or neutri- 
nos, and aware of the uncertainties of the microphysics when considering the interaction between 
radiation and ultradense matter, we introduce a relation between the radiation energy flux den- 
sity and the radiation energy density, i.e. the flux factor, / = T j pr , and the so called variable 
Eddington factor, \ = Pr/ pr, relating the radiation pressure and the radiation energy density. 

There are several of those closure relations reported in the literature (see two recent comprehen- 
sive discussions on this subject in [2TJ Few of them are simply ad hoc relations that smoothly 
interpolate the radiation field between the diffusive and free-streaming regimes. Others, are de- 
rived from a maximum entropy principle [20] or assuming, angular dependence of the radiative 
distribution functions. Even one of them has been motivated from direct transport calculations. 

For our simulations we shall use the Lorentz-Eddington closure relation 

x{f) = l-\^f^P. (33) 

Thus, we can to simulate any radiation phase between diffusion limit and streaming out limit. 
Knowing radiation variables, the field equations ([B])-© together with the anisotropic equation of 
state, the physical variables p, P r , P±, u>, and T can be determined 
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5 Analysis, and some prelimary results 

We integrate the surface equations for 

A(0) = 10.67 Q(0) = -0.008 c(0) = 2.6 (34) 

with 

N = 2.8 £j = 1.0 frj = 0.8 // = 0.001 f u = 1.0 (35) 

and 

M(0) = 2.04 x 10 5 Af o , A.M rac i = O.OlAf (0), i = 15 s and s = 0.25 s (36) 



Radial Pressures at r=C- Radial Pressures at r=C+ 




Figure 2: Hydrodynamic & radiation pressures 

at both sides of the shock front. Notice that both Figure 3: Fluid velocity at r = a(t). Notice that 
pressures have the same order of magnitude and, the boundary surface bounces and expands, 
at the core radiation pressure is larger. 

Our simulations show the energy conditions (specially the strong one) are very sensitive to 
the variation of the local anisotropy. We have found that the possible range for the anisotropic 
parameter is 1 < £ < 4/3, with £ = 1 representing the isotropic limit. Because the energy conditions, 
the behavior of the distribution is different for the same absolute value of the difference of anisotropy 
between the mantle and core |£j — depending on which region is more anisotropic. This was 
found for the case £j = 0.78 and f/j = 0.98 but not for the inverse model ( £j = 0.98 and £// = 0.78). 
In addition, for the same reason a very anisotropic core is not able to maintain a very isotropic 
mantle. 

Another special feature is that anisotropic cores force radiation flux to decay slower than in 
the isotropic case and the velocity of the fluid is slower than of the corresponding isotropic case. 
It seams to indicate that anisotropy slows down the condensation at the core, maintaining it in 
non-ultrarelativistic regimes. This situation has to be further explored in forthcoming studies. 

As it can be appreciated from figure [2j it is found that the hydrodynamic & radiation pressures 
become of same order of magnitude and at the core, radiation pressure is larger. This is consistent 
with the picture we expect for the radiation hydrodynamics within a general relativistic gravitational 
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collapse [35]. Figure ([3]) shows how the shock front bounces, from a collapsing core it turns to expand 
violently. 

The shock strength, N, also appears to be a very significant feature concerning the evolution of 
the distribution, because it drastically increases the fluid velocity behind the shock. From equation 
(j22|) we undertand that this parameter cannot to be less than one because it leads to imaginary 
speed of the shock. The reason is that the fraction When we raise the shock force parameter, the 
fluid velocity increases faster to relativistic values and even for certain critical value (it depends on 
the whole set of initial conditions) it reaches light speed rapidly. 

As we have stressed, these are preliminary results that should deserve further simulations. 
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